#include "CudaEpiModel.h"
#include "CudaParam.h"
#include <thrust/copy.h>
#include <thrust/binary_search.h>
#include <thrust/count.h>
#include <thrust/sort.h>
#include <thrust/remove.h>

bool compareModel(EpiModel& compare_model){
	copyPerson(false);
	for(int i = 0; i < compare_model.nNumPerson; i++){
		if(fabs(compare_model.pvec[i].prs - model->pvec[i].prs) > 1e-5){ 
			cout << "at :" << i <<  " compare: " << compare_model.pvec[i].prs  << " orignial: " << model->pvec[i].prs << endl;
			return false;
		}
	}
	return true;
}

__global__ void tester(){
        for(int i = 0; i < VACCEFFLENGTH + 1; i++){
                printf("%f\n", ddefaultvacceff[i]);
        }
        printf("%d,%d,%d\n",dTARGETCOMMUNITYSIZE,dFAMILIESPERCLUSTER,dWORKGROUPSIZE);
}


void cudaInitialize(EpiModel& input_model){
	model = &input_model;
    // transfer constant  information to constant memory
    cudaMemcpyToSymbol(dvload,model->vload,VLOADNSUB*VLOADNDAY*sizeof(double)); // vload
    cudaMemcpyToSymbol(dseasonality,model->seasonality,MAXRUNLENGTH*sizeof(double)); // seasonality
    cudaMemcpyToSymbol(dfVaccineEfficacyByAge,model->fVaccineEfficacyByAge,TAG*sizeof(double)); // fVaccineEfficacyByAge
    cudaMemcpyToSymbol(dVaccineData,model->VaccineData,NUMVACCINES*sizeof(vaccinedatastruct)); // vaccineData
	cudaMemcpyToSymbol(dcpfc,cpfc,TAG*sizeof(double));
	cudaMemcpyToSymbol(dcpfa,cpfa,TAG*sizeof(double));
	cudaMemcpyToSymbol(dcphcc,cphcc,TAG*sizeof(double));
	cudaMemcpyToSymbol(dcphca,cphca,TAG*sizeof(double));
	cudaMemcpyToSymbol(dwithdrawcdf32,model->withdrawcdf32,3*WITHDRAWDAYS*sizeof(unsigned int)); // withdrawcef32
	cudaMemcpyToSymbol(dincubationcdf32,incubationcdf32,3*sizeof(unsigned int)); // incubationcdf32
	cudaMemcpyToSymbol(dcps,cps,10*sizeof(double)); // cps
	cudaMemcpyToSymbol(d_travel_length_cdf,travel_length_cdf,12*sizeof(double)); // travel_length_cdf
	cudaMemcpyToSymbol(d_travel_pr,travel_pr,TAG*sizeof(double)); // travel_pr
	cudaMemcpyToSymbol(dnSchoolOpeningDays,model->nSchoolOpeningDays,56*sizeof(int));
	cudaMemcpyToSymbol(dnVaccinePriorities2,model->nVaccinePriorities2,PRIORITY_LAST*sizeof(unsigned int));
	
	// copy person vector
	ppl_vec = new dPersonVec(model->nNumPerson);
	copyPerson(true);

	// copy community vector
	comm_vec = new dCommunityVec(model->nNumCommunities);
	copyCommunity(true);

	// copy tract vector
	tract_vec = new dTractVec(model->nNumTractsTotal);
	copyTract(true);

	rand_vec = new dRandVec(model->nNumCommunities);
	//srand(time(NULL));
	srand(1);
	hRandVec rand_temp;
	for(int i = 0; i < model->nNumCommunities; i++){rand_temp.push_back(RandGen(rand(),rand(),rand(),rand()));}
	*rand_vec = rand_temp;

	// initialize community counter and construct comm_woker_vec and comm_visitior_vec
	unsigned int worker_size = 0; 
	unsigned int visitor_size = 0;
	hCommWorkerVec comm_worker_temp,comm_worker_counter_temp;
	hCommVisitorVec comm_visitor_temp,comm_visitor_counter_temp;

	dPerson dummy_ppl;
	for(unsigned int i = 0; i < model->nNumCommunities; i++){
		comm_worker_counter_temp.push_back(dCommWorker(i,0));
		comm_visitor_counter_temp.push_back(dCommVisitor(i,dummy_ppl));

		for(unsigned int j = 0; j < model->commvec[i].workers.size(); j++){
			comm_worker_temp.push_back(dCommWorker(i,model->commvec[i].workers[j]));
		}
		worker_size += model->commvec[i].workers.size();	

		for(list< Person >::iterator it = model->commvec[i].visitors.begin();
		    it != model->commvec[i].visitors.end(); ++it){
			Person &p = *it;
			comm_visitor_temp.push_back(dCommVisitor(i,dPerson(p)));
			visitor_size++;
		}	
	}

	comm_worker_counter = new dCommWorkerVec(model->nNumCommunities);
	comm_visitor_counter = new dCommVisitorVec(model->nNumCommunities);
	*comm_worker_counter = comm_worker_counter_temp;
	*comm_visitor_counter = comm_visitor_counter_temp;

	comm_worker_vec = new dCommWorkerVec(worker_size);
	comm_visitor_vec = new dCommVisitorVec(visitor_size);
	*comm_worker_vec = comm_worker_temp;
	*comm_visitor_vec = comm_visitor_temp;

	// find the start index and end index of worker and visitor for each community
	comm_worker_start_vec = new dIdxVec(model->nNumCommunities);
	comm_worker_end_vec = new dIdxVec(model->nNumCommunities);
	// since visitor is initially empty, so no need to setup
	comm_visitor_start_vec = new dIdxVec(model->nNumCommunities);
	comm_visitor_end_vec = new dIdxVec(model->nNumCommunities);

	
	thrust::lower_bound(
		comm_worker_vec->begin(),comm_worker_vec->end(),
		comm_worker_counter->begin(),comm_worker_counter->end(),
		comm_worker_start_vec->begin());

	thrust::upper_bound(
		comm_worker_vec->begin(),comm_worker_vec->end(),
		comm_worker_counter->begin(),comm_worker_counter->end(),
		comm_worker_end_vec->begin());


	cudaMalloc((void**)&nNumTAPDone,sizeof(unsigned int));
	cudaMalloc((void**)&nNumWantAV,sizeof(unsigned int));
	cudaMalloc((void**)&nNumWantVaccine,sizeof(unsigned int));
	cudaMalloc((void**)&nNumAntiviralsUsed,sizeof(unsigned int));
	
	cudaMemcpy(nNumTAPDone,&model->nNumTAPDone,sizeof(unsigned int),cudaMemcpyHostToDevice);
	cudaMemcpy(nNumWantAV,&model->nNumWantAV,sizeof(unsigned int),cudaMemcpyHostToDevice);
	cudaMemcpy(nNumWantVaccine,&model->nNumWantVaccine,sizeof(unsigned int),cudaMemcpyHostToDevice);
	cudaMemcpy(nNumAntiviralsUsed,&model->nNumAntiviralsUsed,sizeof(unsigned int),cudaMemcpyHostToDevice);
}

void cudaDay(){
	thrust::for_each(ppl_vec->begin(),ppl_vec->end(),
		DayOp1(model->nTimer,model->AVEs,model->AVEi));

	err = cudaThreadSynchronize();
        if(err != cudaSuccess) cout << "cudaThreadSynchronize error: " << cudaGetErrorString(err) << endl;
		
	thrust::for_each(comm_visitor_vec->begin(),comm_visitor_vec->end(),
		DayOp2(model->AVEs,model->AVEi));

	err = cudaThreadSynchronize();
        if(err != cudaSuccess) cout << "cudaThreadSynchronize error: " << cudaGetErrorString(err) << endl;
		
	thrust::device_vector<int> mutex(9,0);
		
	thrust::for_each(comm_vec->begin(),comm_vec->end(),
		DayOp3(thrust::raw_pointer_cast(ppl_vec->data()),
				thrust::raw_pointer_cast(tract_vec->data()),
				thrust::raw_pointer_cast(comm_vec->data()),
				thrust::raw_pointer_cast(rand_vec->data()),
				model->nTimer,model->AVEp,
				model->fSymptomaticAscertainment,model->bTrigger,
				model->fLiberalLeaveCompliance,model->fIsolationCompliance,model->bTravel,
				thrust::raw_pointer_cast(comm_worker_vec->data()),
				thrust::raw_pointer_cast(comm_worker_start_vec->data()),
				thrust::raw_pointer_cast(comm_worker_end_vec->data()),
				thrust::raw_pointer_cast(comm_visitor_vec->data()),
				thrust::raw_pointer_cast(comm_visitor_start_vec->data()),
				thrust::raw_pointer_cast(comm_visitor_end_vec->data()),
				thrust::raw_pointer_cast(mutex.data()))
				);
	
	// DayOp3Improve<<<model->nNumCommunities,500>>>(
		// thrust::raw_pointer_cast(ppl_vec->data()),
		// thrust::raw_pointer_cast(tract_vec->data()),
		// thrust::raw_pointer_cast(comm_vec->data()),
		// thrust::raw_pointer_cast(rand_vec->data()),
		// model->nTimer,model->AVEp,
		// model->fSymptomaticAscertainment,model->bTrigger,
		// model->fLiberalLeaveCompliance,model->fIsolationCompliance,model->bTravel,
		// thrust::raw_pointer_cast(comm_worker_vec->data()),
		// thrust::raw_pointer_cast(comm_worker_start_vec->data()),
		// thrust::raw_pointer_cast(comm_worker_end_vec->data()),
		// thrust::raw_pointer_cast(comm_visitor_vec->data()),
		// thrust::raw_pointer_cast(comm_visitor_start_vec->data()),
		// thrust::raw_pointer_cast(comm_visitor_end_vec->data()),
		// thrust::raw_pointer_cast(mutex.data()));
	
	err = cudaThreadSynchronize();
        if(err != cudaSuccess) cout << "cudaThreadSynchronize error: " << cudaGetErrorString(err) << endl;
}

// true: transfer person vector from host to device
// false: transfer person vector from device to host
void copyPerson(bool host_to_device){
	hPersonVec temp;
	if(host_to_device){

	vector< Person >::iterator pend = model->pvec.end();
	for(vector< Person >::iterator it = model->pvec.begin(); it != pend; it++){
		Person &p = *it;
		dPerson new_ppl(p);
		temp.push_back(new_ppl);
	}
	*ppl_vec = temp;
	
	}else{
	
	temp.resize(model->nNumPerson);
	temp = *ppl_vec;
	for(int i = 0; i < model->nNumPerson; i++){ temp[i].copyToPerson(model->pvec[i]); }
	}
}

// true: transfer community vector from host to device
// false: transfer community vector from device to host
void copyCommunity(bool host_to_device){
	hCommunityVec temp;
	if(host_to_device){

	vector< Community >::iterator cend = model->commvec.end();
	for(vector< Community >::iterator it = model->commvec.begin(); it != cend; it++){
		Community &comm = *it;
		dCommunity new_comm(comm);
		temp.push_back(new_comm);
	}	
	*comm_vec = temp;

	}else{
	
	temp.resize(model->nNumCommunities);
	temp = *comm_vec;
	for(int i = 0; i < model->nNumCommunities; i++){ temp[i].copyToCommunity(model->commvec[i]); }
	}
}

// true: transfer tract vector from host to device
// false: transfer tract vector from device to host
void copyTract(bool host_to_device){
	hTractVec temp;
	if(host_to_device){

	vector< Tract >::iterator cend = model->tractvec.end();
	for(vector< Tract >::iterator it = model->tractvec.begin(); it != cend; it++){
		Tract &t = *it;
		dTract new_t(t);
		temp.push_back(new_t);
	}	
	*tract_vec = temp;

	}else{
	
	temp.resize(model->nNumTractsTotal);
	temp = *tract_vec;
	for(int i = 0; i < model->nNumTractsTotal; i++){ temp[i].copyToTract(model->tractvec[i]); }
	}
}

void cudaTravelEnd(){
	thrust::for_each(comm_visitor_vec->begin(),comm_visitor_vec->end(),
		TravelEndOp(thrust::raw_pointer_cast(ppl_vec->data()),
					thrust::raw_pointer_cast(comm_vec->data())));

	err = cudaThreadSynchronize();
        if(err != cudaSuccess) cout << "cudaThreadSynchronize error: " << cudaGetErrorString(err) << endl;
	
	int size = comm_visitor_vec->size();
	size -= thrust::count_if(comm_visitor_vec->begin(),comm_visitor_vec->end(),TravelRemoveOp());
	thrust::remove_if(comm_visitor_vec->begin(),comm_visitor_vec->end(),TravelRemoveOp());
	comm_visitor_vec->resize(size);
	
	err = cudaThreadSynchronize();
        if(err != cudaSuccess) cout << "cudaThreadSynchronize error: " << cudaGetErrorString(err) << endl;
}

void cudaNight(){
	thrust::device_vector<int> mutex(4,0);
	thrust::device_vector<int> TAP_mutex_1(model->nNumCommunities,0);
	thrust::device_vector<int> TAP_mutex_2(1,0);
	
	thrust::for_each(comm_vec->begin(),comm_vec->end(),
		NightOp(thrust::raw_pointer_cast(ppl_vec->data()),
				thrust::raw_pointer_cast(tract_vec->data()),
				thrust::raw_pointer_cast(comm_vec->data()),
				thrust::raw_pointer_cast(rand_vec->data()),
				
				model->nTimer,model->AVEp,
				model->fSymptomaticAscertainment,model->bTrigger,
				model->fLiberalLeaveCompliance,model->fIsolationCompliance,
				model->fQuarantineCompliance,model->bTravel,model->nAscertainmentDelay,
				model->fContactAscertainment,
				
				nNumTAPDone,nNumWantAV,
				
				thrust::raw_pointer_cast(comm_visitor_vec->data()),
				thrust::raw_pointer_cast(comm_visitor_start_vec->data()),
				thrust::raw_pointer_cast(comm_visitor_end_vec->data()),
				
				thrust::raw_pointer_cast(mutex.data()),
				thrust::raw_pointer_cast(TAP_mutex_1.data()),
				thrust::raw_pointer_cast(TAP_mutex_2.data()))
				);
				
	err = cudaThreadSynchronize();
        if(err != cudaSuccess) cout << "cudaThreadSynchronize error: " << cudaGetErrorString(err) << endl;
		
}

void cudaTravelStart(){
	thrust::device_vector<dCommVisitor> new_visitor_vec(model->nNumPerson);

	thrust::for_each(ppl_vec->begin(),ppl_vec->end(),
		TravelStartOp(thrust::raw_pointer_cast(ppl_vec->data()),
					  thrust::raw_pointer_cast(tract_vec->data()),
					  thrust::raw_pointer_cast(comm_vec->data()),
					  thrust::raw_pointer_cast(rand_vec->data()),
					  thrust::raw_pointer_cast(new_visitor_vec.data()),
					  model->nNumTractsTotal));
	
	int size = thrust::count_if(new_visitor_vec.begin(),new_visitor_vec.end(),TravelRemoveOp());
	size = model->nNumPerson - size;
	thrust::remove_if(new_visitor_vec.begin(),new_visitor_vec.end(),TravelRemoveOp());
	new_visitor_vec.resize(size);
	
	int comm_visitor_size = comm_visitor_vec->size();
	int new_comm_visitor_size = comm_visitor_size + size;
	comm_visitor_vec->resize(new_comm_visitor_size);
	
	thrust::copy(new_visitor_vec.begin(),new_visitor_vec.end(),comm_visitor_vec->begin()+comm_visitor_size);
	
	thrust::sort(comm_visitor_vec->begin(),comm_visitor_vec->end());
	
	thrust::lower_bound(
		comm_visitor_vec->begin(),comm_visitor_vec->end(),
		comm_visitor_counter->begin(),comm_visitor_counter->end(),
		comm_visitor_start_vec->begin());

	thrust::upper_bound(
		comm_visitor_vec->begin(),comm_visitor_vec->end(),
		comm_visitor_counter->begin(),comm_visitor_counter->end(),
		comm_visitor_end_vec->begin());
}

void cudaResponse(){
	thrust::for_each(tract_vec->begin(),tract_vec->end(),ActOp0(model->nTimer));
	
	err = cudaThreadSynchronize();
        if(err != cudaSuccess) cout << "cudaThreadSynchronize error: " << cudaGetErrorString(err) << endl;

	// 1. Count cumulative number of ascertained cases
	if(!model->bTrigger){
		thrust::device_vector<int> comm_nNumAscertained_vec(model->nNumCommunities);
		thrust::transform(comm_vec->begin(),comm_vec->end(),comm_nNumAscertained_vec.begin(),ActOp1());
		int nNumAscertained = thrust::reduce(comm_nNumAscertained_vec.begin(),comm_nNumAscertained_vec.end());
	
		if (nNumAscertained/(double)model->nNumPerson>model->fResponseThreshhold) {
			model->bTrigger=true;
			model->nTriggerTime=model->nTimer+model->nTriggerDelay*2;
			cout << "Response starts on day " << (model->nTriggerTime/2) << ": ascertained: " 
			<< nNumAscertained << "/" << model->nNumPerson << "=" << nNumAscertained/(double)model->nNumPerson 
			<< " on day " << (model->nTimer/2) << endl;
		}
	}
	
	err = cudaThreadSynchronize();
        if(err != cudaSuccess) cout << "cudaThreadSynchronize error: " << cudaGetErrorString(err) << endl;
	
	// 2. Epidemic is detected, start reactive vaccination strategies
	if(model->bTrigger && model->nTimer==model->nTriggerTime) {

		thrust::for_each(tract_vec->begin(),tract_vec->end(),
			ActOp2(thrust::raw_pointer_cast(ppl_vec->data()),
				thrust::raw_pointer_cast(comm_vec->data()),
				nNumWantVaccine,
				model->eAVPolicy,model->fQuarantineCompliance,model->nSchoolClosureDays,
				model->nSchoolClosurePolicy,model->eVaccinationStrategy,model->nAscertainmentDelay));
			
	}
	
	err = cudaThreadSynchronize();
        if(err != cudaSuccess) cout << "cudaThreadSynchronize error: " << cudaGetErrorString(err) << endl;
	
	 // 3. change vaccine priorities if needed
	if (model->nPriorityChangeTime>0 && model->nPriorityChangeTime==model->nTimer/2) {
		thrust::for_each(ppl_vec->begin(),ppl_vec->end(),
			ActOp3(thrust::raw_pointer_cast(ppl_vec->data()),
				   model->nNumPerson,model->nHighestPriority
				   ));
		
	}
	
	err = cudaThreadSynchronize();
        if(err != cudaSuccess) cout << "cudaThreadSynchronize error: " << cudaGetErrorString(err) << endl;
	
	// 4. Ring vaccination and consider closing or re-opening schools
	if (model->bTrigger && model->nTimer>=model->nTriggerTime) {
		if (model->eVaccinationStrategy==RINGVACTRACT) {
	
			thrust::for_each(tract_vec->begin(),tract_vec->end(),
				ActOp41(thrust::raw_pointer_cast(ppl_vec->data()),
						thrust::raw_pointer_cast(comm_vec->data()),
						nNumWantVaccine,
						model->nAscertainmentDelay));
						
		}
		
		if (model->nSchoolClosureDays>0) {
			if (model->nSchoolClosurePolicy==2) { // close schools by ascertainment
				thrust::for_each(tract_vec->begin(),tract_vec->end(),
					ActOp42(thrust::raw_pointer_cast(ppl_vec->data()),
							thrust::raw_pointer_cast(comm_vec->data()),
							model->nAscertainmentDelay));
				// original program has output string
			}
			
			thrust::for_each(tract_vec->begin(),tract_vec->end(),ActOp43(model->nTimer));
			// original program has output string
		}
	}
	
	err = cudaThreadSynchronize();
        if(err != cudaSuccess) cout << "cudaThreadSynchronize error: " << cudaGetErrorString(err) << endl;
	
	// 5. Count available vaccines
	for (int i=0; i<model->nNumVaccineTypes; i++)
		if (model->vaccineproductionschedule[i][model->nTimer/2]>0)
			model->nVaccineSupply[i] += model->vaccineproductionschedule[i][model->nTimer/2];  // new vaccine arrives

	unsigned int nVacLim[NUMVACCINES];
	unsigned int totalvaccinesupply=0;
	
	for (int i=0; i<model->nNumVaccineTypes; i++) {
		nVacLim[i] = model->nVaccineSupply[i]-model->nNumVaccineDosesUsed[i]; // how many vaccine doses are available today?
		if (model->nVaccineSupply[i]<=model->nNumVaccineDosesUsed[i]) // because we are using unsigned ints
			nVacLim[i]=0;
		totalvaccinesupply += nVacLim[i];
	}
	if (totalvaccinesupply>model->nVaccineDailyLimit) { // limited overall vaccine daily delivery capacity
		double factor = model->nVaccineDailyLimit/(double)totalvaccinesupply;
		for (int i=0; i<model->nNumVaccineTypes; i++) nVacLim[i] *= factor;
	}
	
	err = cudaThreadSynchronize();
        if(err != cudaSuccess) cout << "cudaThreadSynchronize error: " << cudaGetErrorString(err) << endl;
	
	// 6. Distribute vaccines
	thrust::device_vector<unsigned int> nVacLim_vec(model->nNumPerson);
	thrust::device_vector<unsigned int> d_nNumVaccineDosesUsed(
		model->nNumVaccineDosesUsed,model->nNumVaccineDosesUsed+NUMVACCINES);
	
	
	for (int vacnum=0; vacnum<model->nNumVaccineTypes; vacnum++) {
		
		if (nVacLim[vacnum]>0) {
			double vFrac=0.0;
			
			thrust::device_vector<unsigned int> wantthisvac(PRIORITY_LAST,0);
			thrust::for_each(ppl_vec->begin(),ppl_vec->end(),
				ActOp61(thrust::raw_pointer_cast(wantthisvac.data()),vacnum));
				
			for (int priority=1; priority<=model->nHighestPriority; priority++) {
			
				unsigned int wantvac = wantthisvac[priority];
				if (wantvac > 0) {
					vFrac = nVacLim[vacnum]/(double)wantvac;
				}
				thrust::device_vector<int> mutex(1,0);
				thrust::for_each(ppl_vec->begin(),ppl_vec->end(),
					ActOp62(priority,vacnum,vFrac,
							thrust::raw_pointer_cast(rand_vec->data()),
							thrust::raw_pointer_cast(d_nNumVaccineDosesUsed.data()),
							thrust::raw_pointer_cast(nVacLim_vec.data()),
							nNumWantVaccine,
							model->ePrevaccinationStrategy,
							thrust::raw_pointer_cast(mutex.data())
							));
				
				unsigned int nVacLim_sum = thrust::reduce(nVacLim_vec.begin(),nVacLim_vec.end());
				if(nVacLim[vacnum] > nVacLim_sum) nVacLim[vacnum] -= nVacLim_sum;
				else nVacLim[vacnum] = 0;
				
				err = cudaThreadSynchronize();
					if(err != cudaSuccess) cout << "cudaThreadSynchronize error: " << cudaGetErrorString(err) << endl;
			}
		}
	}
	
	err = cudaThreadSynchronize();
        if(err != cudaSuccess) cout << "cudaThreadSynchronize error: " << cudaGetErrorString(err) << endl;
	
	// 7. Count available antiviral courses
	cudaMemcpy(&model->nNumAntiviralsUsed,nNumAntiviralsUsed,sizeof(unsigned int),cudaMemcpyHostToDevice);
	cudaMemcpy(&model->nNumWantAV,nNumWantAV,sizeof(unsigned int),cudaMemcpyDeviceToHost);
	
	unsigned int uAVFraction;
	unsigned int nAVLim = min(model->nAVTotalLimit-model->nNumAntiviralsUsed, model->nAVDailyLimit);
	if (model->nAVTotalLimit <= model->nNumAntiviralsUsed) nAVLim=0;
	
	
	if (model->nNumWantAV <= nAVLim)
		uAVFraction = UINT_MAX;
	else
		uAVFraction = (unsigned int)(nAVLim/(double)model->nNumWantAV*UINT_MAX);
	
	thrust::for_each(ppl_vec->begin(),ppl_vec->end(),
		ActOp71(nNumAntiviralsUsed,nNumWantAV,uAVFraction,
			thrust::raw_pointer_cast(rand_vec->data())));
	
	cudaMemcpy(&model->nNumAntiviralsUsed,nNumAntiviralsUsed,sizeof(unsigned int),cudaMemcpyDeviceToHost);
	
	for (int i=0; i<model->nNumVaccineTypes; i++)
		model->nNumVaccineDosesUsedHistory[i].push_back(d_nNumVaccineDosesUsed[i]);
	model->nNumAntiviralsUsedHistory.push_back(model->nNumAntiviralsUsed);
	
	err = cudaThreadSynchronize();
        if(err != cudaSuccess) cout << "cudaThreadSynchronize error: " << cudaGetErrorString(err) << endl;
	
}

void cudaStats(){
	unsigned int nNumSymptomatic=0;
    unsigned int nNumCumulativeSymptomatic=0;
	thrust::device_vector<unsigned int> results(model->nNumCommunities);
	
	thrust::transform(comm_vec->begin(),comm_vec->end(),results.begin(),StatsOp1());
	nNumSymptomatic = reduce(results.begin(),results.end());
	
	err = cudaThreadSynchronize();
        if(err != cudaSuccess) cout << "cudaThreadSynchronize error: " << cudaGetErrorString(err) << endl;
	
	thrust::transform(comm_vec->begin(),comm_vec->end(),results.begin(),StatsOp2());
	nNumCumulativeSymptomatic = reduce(results.begin(),results.end());
	
	err = cudaThreadSynchronize();
        if(err != cudaSuccess) cout << "cudaThreadSynchronize error: " << cudaGetErrorString(err) << endl;
	
	model->nNumCumulativeSymptomaticHistory.push_back(nNumCumulativeSymptomatic);
    model->nNumSymptomaticHistory.push_back(nNumSymptomatic);
	
	cout << "nNumCumulativeSymptomatic: " << nNumCumulativeSymptomatic << endl;
	cout << "nNumSymptomatic: " << nNumSymptomatic << endl;
	
}

void clean(){
	delete ppl_vec;
	delete comm_vec;
	delete comm_worker_vec; 
	delete comm_worker_start_vec; 
	delete comm_worker_end_vec;  

	delete comm_visitor_vec; 
	delete comm_visitor_start_vec; 
	delete comm_visitor_end_vec; 

	delete comm_worker_counter;
	delete comm_visitor_counter;
	delete rand_vec;
}


__device__ void nightInfectSusceptibles(
	const dPerson& infected,dCommunity& comm,
	dPerson* ppl_vec,dCommunity* comm_vec,RandGen* rand_vec,
	int timer,	
	double AVEp,
	double fSymptomaticAscertainment,
	int bTrigger,
	double fLiberalLeaveCompliance,
	double fIsolationCompliance,
	bool bTravel,
	dCommVisitor* comm_visitor_vec, unsigned int comm_visitor_start_idx, unsigned int comm_visitor_end_idx,
	int* mutex, int target_value){
	
	const double *cpf = (dPerson::isChild(infected)?dcpfc:dcpfa); // probability of family transmission
	const double *cphc = (dPerson::isChild(infected)?dcphcc:dcphca); // probability of household cluster transmission
	
	for (unsigned int pid2=comm.nFirstPerson; pid2<comm.nLastPerson; pid2++) {
		dPerson &p2 = ppl_vec[pid2];
		if (dPerson::isSusceptible(p2) && p2.nTravelTimer<=0) { // && p.id!=p2.id) {
			if (infected.family==p2.family){
				InfectOp new_infect1(p2,infected,cpf[p2.age],FROMFAMILYNIGHT,
									rand_vec[comm.id],timer,AVEp,
									fSymptomaticAscertainment,bTrigger,fLiberalLeaveCompliance,
									fIsolationCompliance,comm_vec);
			
				if (new_infect1.infect()) continue;
			}
			if (!dPerson::isWithdrawn(infected) && 
			!dPerson::isQuarantined(infected) && !dPerson::isQuarantined(p2)) {
				InfectOp new_infect2(p2,infected,comm.cpcm[p2.age],FROMCOMMUNITYNIGHT,
									rand_vec[comm.id],timer,AVEp,
									fSymptomaticAscertainment,bTrigger,fLiberalLeaveCompliance,
									fIsolationCompliance,comm_vec);
									
				if (!new_infect2.infect()) { // transmission within home community
					if (infected.nHomeNeighborhood==p2.nHomeNeighborhood){ // transmission within neighborhood
						InfectOp new_infect3(p2,infected,comm.cpnh[p2.age],FROMNEIGHBORHOODNIGHT,
											rand_vec[comm.id],timer,AVEp,
											fSymptomaticAscertainment,bTrigger,fLiberalLeaveCompliance,
											fIsolationCompliance,comm_vec);
					
						if (new_infect3.infect()) continue;
					}
					if (infected.householdcluster==p2.householdcluster){ // transmission within household cluster
						InfectOp new_infect4(p2,infected,cphc[p2.age],FROMHHCLUSTERNIGHT,
											rand_vec[comm.id],timer,AVEp,
											fSymptomaticAscertainment,bTrigger,fLiberalLeaveCompliance,
											fIsolationCompliance,comm_vec);
						new_infect4.infect();
					}
				}
			}
		}
	}
	
	Barrier(mutex,target_value);
	
	// check for susceptible travelers
	if (bTravel) {
  
		for (unsigned int i = comm_visitor_start_idx; i < comm_visitor_end_idx; i++) {
			dPerson &p2 = comm_visitor_vec[i].visitor;
			if (dPerson::isSusceptible(p2)) {
				if (infected.family==p2.family){  // transmission within household
					InfectOp new_infect1(p2,infected,cpf[p2.age],FROMFAMILY,
										rand_vec[comm.id],timer,AVEp,
										fSymptomaticAscertainment,bTrigger,fLiberalLeaveCompliance,
										fIsolationCompliance,comm_vec);
					
					if (new_infect1.infect()) continue;
				}
				if (!dPerson::isWithdrawn(infected) 
				&& !dPerson::isQuarantined(infected) && !dPerson::isQuarantined(p2)) {
					InfectOp new_infect2(p2,infected,comm.cpcm[p2.age],FROMCOMMUNITYNIGHT,
										rand_vec[comm.id],timer,AVEp,
										fSymptomaticAscertainment,bTrigger,fLiberalLeaveCompliance,
										fIsolationCompliance,comm_vec);
					
					if (!new_infect2.infect()) { // transmission within home community
						if (infected.nHomeNeighborhood==p2.nHomeNeighborhood){ // transmission within neighborhood
							InfectOp new_infect3(p2,infected,comm.cpnh[p2.age],FROMNEIGHBORHOODNIGHT,
												rand_vec[comm.id],timer,AVEp,
												fSymptomaticAscertainment,bTrigger,fLiberalLeaveCompliance,
												fIsolationCompliance,comm_vec);

							if (new_infect3.infect()) continue;
						}
						if (infected.householdcluster==p2.householdcluster){ // transmission within household cluster
							InfectOp new_infect4(p2,infected,cphc[p2.age],FROMHHCLUSTERNIGHT,
												rand_vec[comm.id],timer,AVEp,
												fSymptomaticAscertainment,bTrigger,fLiberalLeaveCompliance,
												fIsolationCompliance,comm_vec);
							
							new_infect4.infect();
						}
					}
				}
			}
		}
	}
	
	Barrier(mutex+1,target_value);
}
